clusteredSE<-function(model,data,cluster){
  nona.x<-na.omit(data[,c(cluster, all.vars(formula(model)))])
  M <- length(unique(data[[cluster]])) # Number of groups   
  N <- length(nona.x[[cluster]]) # Number of observations
  if(class(model)=="zeroinfl"){K <- dim(model$vcov)[1]} else {K <- model$rank} # Number of parameters                 
  dfc <- (M/(M-1))*((N-1)/(N-K)) #
  uj  <- apply(estfun(model),2,function(x) tapply(x,nona.x[[cluster]], sum))
  if(class(model)=="zeroinfl"){
    vcovCL <- dfc[1]*sandwich(model, meat=crossprod(uj)/N)*1
    list(summary=coeftest(model, vcovCL))$summary} 
  else {
    vcovCL<-dfc*sandwich(model, meat=crossprod(uj)/N)
    clus.se<-coeftest(model,vcovCL)               
    return(list(clus.se))}}

coef_zinb<-function(model,label='',remove.factors=T,drop=''){
  coefs<-tibble(term=model@coef.names,
                estimate=as.numeric(model@coef[1:length(model@coef.names)]),
                std.error=as.numeric(model@se[1:length(model@coef.names)[1]]),
                model=label)
  coefs$term<-str_remove(coefs$term,'count_') # rename coefficient of count equation
  coefs<-subset(coefs,str_count(coefs$term,'zero_')==0) # drop coefficient of zero-inflation equation
  if (remove.factors==T){coefs<-subset(coefs,str_count(coefs$term,'factor\\(')==0)}
  if (drop!=''){for(i in 1:length(drop)){coefs<-subset(coefs,coefs$term!=drop[i])}}
  return(coefs)}

corstarsl <- function(x){
  require(Hmisc)
  x <- as.matrix(x)
  R <- rcorr(x)$r
  p <- rcorr(x)$P
  mystars <- ifelse(p < .01, "***", ifelse(p < .05, "** ", ifelse(p < .1, "* ", " "))) ## define notions for significance levels; spacing is important.
  R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1] ## trunctuate the matrix that holds the correlations to two decimal
  ## build a new matrix that includes the correlations with their apropriate stars
  Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x))
  diag(Rnew) <- paste(diag(R), " ", sep="")
  rownames(Rnew) <- colnames(x)
  colnames(Rnew) <- paste(colnames(x), "", sep="")
  ## remove upper triangle
  Rnew <- as.matrix(Rnew)
  Rnew[upper.tri(Rnew, diag = TRUE)] <- ""
  Rnew <- as.data.frame(Rnew)
  Rnew <- cbind(Rnew[1:length(Rnew)-1]) ## remove last column and return the matrix (which is now a data frame)
  return(Rnew)}